# Load required packages
library(cjoint)

# Load the data
load("survey_data_nominees_final.RData")

#########################
# Figure 1: Main result #
#########################
# Subsetting to respondents who answered the support question
support.dataset <- survey.dataset[which(!is.na(survey.dataset$support_nominee)),]

# Default baselines
baselines_list <- list()

# Estimate the AMCEs
results.support <- amce(support_nominee ~ judge_age + judge_race + judge_law_school + judge_job + judge_politics
                        + judge_transgender + judge_sexuality + judge_rhetoric, data=support.dataset, 
                        baselines = baselines_list, cluster=T, respondent.id = "ResponseId")
summary(results.support)

# Estimate the AMCEs, results by respondent partisanship
results.support.by.party <- amce(support_nominee ~ judge_age*Party + judge_race*Party 
                                 + judge_law_school*Party + judge_job*Party 
                                 + judge_politics*Party + judge_transgender*Party 
                                 + judge_sexuality*Party+ judge_rhetoric*Party,
                                 data=support.dataset, respondent.varying="Party", baselines = baselines_list, cluster=T, 
                                 respondent.id = "ResponseId")
summary(results.support.by.party)

### Figure 1: Plot of differences when moving from baseline ###
# Values come from above models
# All respondents first, then D, then R
coefs_from_baseline_sexuality <- c(-0.08096516,
                                   -0.0562235,
                                   -0.0979756)
# All respondents first, then D, then R
coefs_from_baseline_gender <- c(-0.13838154,
                                -0.1225160,
                                -0.1544877)

# Confidence intervals around the differences
sexuality_se <- c(0.019370,0.025061,0.024170)
gender_se <- c(0.020081,0.025880,0.023581)

sexuality_ci_low <- coefs_from_baseline_sexuality - (sexuality_se*1.96)
sexuality_ci_high <- coefs_from_baseline_sexuality + (sexuality_se*1.96)
sexuality_ci_low_90 <- coefs_from_baseline_sexuality - (sexuality_se*1.65)
sexuality_ci_high_90 <- coefs_from_baseline_sexuality + (sexuality_se*1.65)

gender_ci_low <- coefs_from_baseline_gender - (gender_se*1.96)
gender_ci_high <- coefs_from_baseline_gender + (gender_se*1.96)
gender_ci_low_90 <- coefs_from_baseline_gender - (gender_se*1.65)
gender_ci_high_90 <- coefs_from_baseline_gender + (gender_se*1.65)

# Plotting
setEPS(width=4,height=4)
postscript("figures/figure1-final.eps")
par(mar=c(4.5,4.5,1,1), mgp=c(2,.25,0))
plot(coefs_from_baseline_sexuality,c(0.5,0.35,0.2),ylim=c(0,2),xlim=c(-.3,.1), tck=0, cex.main=1.4, main="", cex.lab=0.75, cex.axis=0.75, 
     pch=c(16,17,15), cex=2, xlab="Change in Support from Baseline (Cisgender or Straight)", xaxt="n", yaxt="n", ylab="", lwd=1.5, col=c("black","blue","red"))
abline(v=0, lty=2)
axis(1,at=seq(-.3,.1,by=.1), cex.axis=0.75, 
     labels=c("-.3","-.2","-.1","0",".1"), tck=0.01, las=1)
axis(2, at=c(0.35,1.35),cex.axis=0.75,labels=c("Gay/Lesbian","Transgender"), padj=0, las=1, tick=F)
legend("topleft", c("All Respondents","Democrats","Republicans"), col=c("black","blue","red"), pch=c(16,17,15), cex=.7)
segments(sexuality_ci_low, c(0.5,0.35,0.2), sexuality_ci_high, c(0.5,0.35,0.2), lwd=1, col=c("black","blue","red"))
segments(sexuality_ci_low_90, c(0.5,0.35,0.2), sexuality_ci_high_90, c(0.5,0.35,0.2), lwd=3, col=c("black","blue","red"))

points(x=coefs_from_baseline_gender, y=c(1.5,1.35,1.2), pch=c(16,17,15), col=c("black","blue","red"), cex=2)
segments(gender_ci_low, c(1.5,1.35,1.2), gender_ci_high, c(1.5,1.35,1.2), lwd=1, col=c("black","blue","red"))
segments(gender_ci_low_90, c(1.5,1.35,1.2), gender_ci_high_90, c(1.5,1.35,1.2), lwd=3, col=c("black","blue","red"))
dev.off()

###############################################################
# In-text discussion of coefficients associated with Figure 1 #
###############################################################
### When averaging across all respondents, prospective
### nominees who identify as gay or lesbian experience an 8.1 percentage point (p<0.001) decrease
### in support compared to straight nominees, and transgender nominees face a 13.8 percentage
### point (p<0.001) decrease in comparison to cisgender nominees. 

# Gay/lesbian penalty, coefficient and p-value
round(summary(results.support)$amce[14,3],3) # coefficient
round(summary(results.support)$amce[14,6],3) # p-value
# Transgender penalty, coefficient and p-value
round(summary(results.support)$amce[15,3],3) # coefficient
round(summary(results.support)$amce[15,6],3) # p-value

### Judges who graduated from non-top 100 schools received 12.3 percentage points (p$<$0.001) 
### lower support in comparison to elite Ivy graduates.

# Law school outside top 100, coefficient and p-value
round(summary(results.support)$amce[7,3],3) # coefficient
round(summary(results.support)$amce[7,6],3) # p-value

### Similar to previous studies, the negative effect observed for transgender 
### individuals is significantly larger than that of gays and lesbians (p=0.037).

# Test of difference in gay and transgender coefficients
z_gt_difference <- (results.support$estimates$judgesexuality[1]-results.support$estimates$judgetransgender[1])/
  (sqrt(((results.support$estimates$judgesexuality[2])^2)+((results.support$estimates$judgetransgender[2])^2)-2*results.support$vcov.prof["judgetransgenderTransgender","judgesexualityGayorlesbian"]))
round(2*pnorm(z_gt_difference, lower=F),3)

### Although baseline support for Biden’s prospective nominee is considerably higher
### among Democrats than Republicans (by 53.5 p.p.) 
round(mean(support.dataset$support_nominee[which(support.dataset$Party == "democrat")]) - mean(support.dataset$support_nominee[which(support.dataset$Party == "republican")]),3)

### For instance, the decrease in support for transgender nominees is 15.4 percentage points
### (p<0.001) among Republicans and 12.3 percentage points (p<0.001) among Democrats. 
### These treatment effects are not distinguishable from one another (p=0.361). 

# R coefficient and p-value
round(summary(results.support.by.party)$Party3amce[13,3],3) # coefficient
round(summary(results.support.by.party)$Party3amce[13,6],3) # p-value
# D coefficient and p-value
round(summary(results.support.by.party)$Party1amce[13,3],3) # coefficient
round(summary(results.support.by.party)$Party1amce[13,6],3) # p-value

# Calculating difference in R and D (values from the cjoint output)
# Transgender
r.d.transgender.difference <- results.support.by.party$cond.estimates$`judgetransgender:Party`[1,2]
r.d.transgender.se <- results.support.by.party$cond.estimates$`judgetransgender:Party`[2,2]
r.d.transgender.z.value <- abs(r.d.transgender.difference/r.d.transgender.se)
round(2*pnorm(r.d.transgender.z.value, lower=F),3)

### Republican respondents favor moderate nominees over very liberal nominees by 14.0 percentage
### points (p<0.001), the effects of ideology for Democrats are indistinguishable from zero at p=0.050.

# Moderate ideology, R coefficient and p-value
round(summary(results.support.by.party)$Party3amce[12,3],3) # coefficient
round(summary(results.support.by.party)$Party3amce[12,6],3) # p-value

# Moderate ideology, D coefficient and p-value
round(summary(results.support.by.party)$Party1amce[12,3],3) # coefficient
round(summary(results.support.by.party)$Party1amce[12,6],3) # p-value
# Somewhat liberal ideology, D coefficient and p-value
round(summary(results.support.by.party)$Party1amce[11,3],3) # coefficient
round(summary(results.support.by.party)$Party1amce[11,6],3) # p-value
# Liberal ideology, D coefficient and p-value
round(summary(results.support.by.party)$Party1amce[10,3],3) # coefficient
round(summary(results.support.by.party)$Party1amce[10,6],3) # p-value

#############################
# Figure 2a: Double penalty #
#############################
# Creating a variable for both/one/neither transgender and gay treatment
support.dataset$judge_combined_transgender_gay <- ifelse(support.dataset$judge_transgender == "Transgender" & support.dataset$judge_sexuality == "Gay or lesbian", "Transgender and gay",
                                                               ifelse((support.dataset$judge_transgender == "Transgender" & support.dataset$judge_sexuality == "Straight") | 
                                                                        (support.dataset$judge_transgender == "Not transgender" & support.dataset$judge_sexuality == "Gay or lesbian"), "Just transgender or just gay","Cisgender and straight"))
support.dataset$judge_combined_transgender_gay <- as.factor(support.dataset$judge_combined_transgender_gay)

# Setting baselines
baselines_list <- list()
# Making just transgender or just gay the baseline
baselines_list$judge_combined_transgender_gay <- "Just transgender or just gay"

# Estimate the AMCEs for gay/transgender double penalty 
results.both.trans.gay.support.pooled <- amce(support_nominee ~ judge_age + judge_race 
                                              + judge_law_school + judge_job 
                                              + judge_politics + judge_combined_transgender_gay
                                              + judge_rhetoric,
                                              data=support.dataset, baselines = baselines_list, cluster=T, 
                                              respondent.id = "ResponseId")
summary(results.both.trans.gay.support.pooled)

### Plot of differences when moving from baseline ###
# Values come from above model
coef_from_baseline_double_pooled <- c(-0.0604171)
coef_from_baseline_neither_pooled <- c(0.1580701)

# Confidence intervals around the differences
double_se_pooled <- c(0.023694)
neither_se_pooled <- c(0.024604)

double_ci_low_pooled <- coef_from_baseline_double_pooled - (double_se_pooled*1.96)
double_ci_high_pooled <- coef_from_baseline_double_pooled + (double_se_pooled*1.96)
double_ci_low_pooled_90 <- coef_from_baseline_double_pooled - (double_se_pooled*1.65)
double_ci_high_pooled_90 <- coef_from_baseline_double_pooled + (double_se_pooled*1.65)

neither_ci_low_pooled <- coef_from_baseline_neither_pooled - (neither_se_pooled*1.96)
neither_ci_high_pooled <- coef_from_baseline_neither_pooled + (neither_se_pooled*1.96)
neither_ci_low_pooled_90 <- coef_from_baseline_neither_pooled - (neither_se_pooled*1.65)
neither_ci_high_pooled_90 <- coef_from_baseline_neither_pooled + (neither_se_pooled*1.65)

# Plotting
setEPS(width=4,height=4)
postscript("figures/figure2a-final.eps")
par(mar=c(4.5,4.5,1,1), mgp=c(2,.25,0))
plot(coef_from_baseline_double_pooled,c(0.35),ylim=c(0,2),xlim=c(-.2,.3), tck=0, cex.main=1.4, main="", cex.lab=0.7, cex.axis=0.75, 
     pch=c(16), cex=2, xlab="Change in Support from Baseline (Just Gay or Just Transgender)", xaxt="n", yaxt="n", ylab="", lwd=1.5, col=c("black"))
abline(v=0, lty=2)
axis(1,at=seq(-.2,.3,by=.1), cex.axis=0.75, 
     labels=c("-.2","-.1","0",".1",".2",".3"), tck=0.01, las=1)
axis(2, at=c(0.3,1.3),cex.axis=0.75,labels=c("Gay and 
Transgender","Straight and 
                                              Cisgender"), padj=0, las=1, tick=F)
legend("topleft", c("All Respondents"), col=c("black"), pch=c(16), cex=.7)
segments(double_ci_low_pooled, c(0.35), double_ci_high_pooled, c(0.35), lwd=1, col=c("black"))
segments(double_ci_low_pooled_90, c(0.35), double_ci_high_pooled_90, c(0.35), lwd=3, col=c("black"))

points(x=coef_from_baseline_neither_pooled, y=c(1.35), pch=c(16), col=c("black"), cex=2)
segments(neither_ci_low_pooled, c(1.35), neither_ci_high_pooled, c(1.35), lwd=1, col=c("black"))
segments(neither_ci_low_pooled_90, c(1.35), neither_ci_high_pooled_90, c(1.35), lwd=3, col=c("black"))
dev.off()

################################################################
# In-text discussion of coefficients associated with Figure 2a #
################################################################
### In Figure 2a, we illustrate that judges who are both transgender and gay pay a greater
### penalty (6.0 percentage point lower support, p=0.011) than judges with just one of those traits. 
# Gay/transgender double penalty, coefficient and p-value
round(summary(results.both.trans.gay.support.pooled)$amce[4,3],3) # coefficient
round(summary(results.both.trans.gay.support.pooled)$amce[4,6],3) # p-value

########################################
# Figure 2b: Greater penalty for women #
########################################
# Subsetting to just men or transgender men nominees
support.dataset.just.men <- support.dataset[which(support.dataset$judge_gender %in% c("Man","Transgender man")),]
# Revised gender variable
support.dataset.just.men$judge_gender_revised <- as.factor(ifelse(support.dataset.just.men$judge_gender == "Man","Man","Transgender man"))

# Setting baselines
baselines_list <- list()

# Estimate the AMCEs for all respondents, just men or transgender men nominees
results.support.just.men <- amce(support_nominee ~ judge_age + judge_race + judge_law_school + judge_job + judge_politics
                                 + judge_gender_revised + judge_sexuality + judge_rhetoric, data=support.dataset.just.men, 
                                 baselines = baselines_list, cluster=T, respondent.id = "ResponseId")
summary(results.support.just.men)

# Subsetting to just women or transgender women nominees
support.dataset.just.women <- support.dataset[which(support.dataset$judge_gender %in% c("Woman","Transgender woman")),]
# Revised gender variable
support.dataset.just.women$judge_gender_revised <- as.factor(ifelse(support.dataset.just.women$judge_gender == "Woman","Woman","Transgender woman"))
support.dataset.just.women$judge_gender_revised <- relevel(support.dataset.just.women$judge_gender_revised, ref="Woman")

# Setting baselines
baselines_list <- list()

# Estimate the AMCEs for all respondents, just women or transgender women nominees
results.support.just.women <- amce(support_nominee ~ judge_age + judge_race + judge_law_school + judge_job + judge_politics
                                   + judge_gender_revised + judge_sexuality + judge_rhetoric, data=support.dataset.just.women, 
                                   baselines = baselines_list, cluster=T, respondent.id = "ResponseId")
summary(results.support.just.women)

### Plot of differences when moving from baseline ###
# Values come from above models, gay/lesbian, man first then woman
coef_from_baseline_gay_mw <- c(-0.0334528,-0.1196301)
# Transgender, man first then woman
coef_from_baseline_trans_mw <- c(-0.1039549,-0.1681376)

# Confidence intervals around the differences
gay_mw_se <- c(0.028190,0.026813)
trans_mw_se <- c(0.028724,0.027809)

gay_mw_ci_low_pooled <- coef_from_baseline_gay_mw - (gay_mw_se*1.96)
gay_mw_ci_high_pooled <- coef_from_baseline_gay_mw + (gay_mw_se*1.96)
gay_mw_ci_low_pooled_90 <- coef_from_baseline_gay_mw - (gay_mw_se*1.65)
gay_mw_ci_high_pooled_90 <- coef_from_baseline_gay_mw + (gay_mw_se*1.65)

trans_mw_ci_low_pooled <- coef_from_baseline_trans_mw - (trans_mw_se*1.96)
trans_mw_ci_high_pooled <- coef_from_baseline_trans_mw + (trans_mw_se*1.96)
trans_mw_ci_low_pooled_90 <- coef_from_baseline_trans_mw - (trans_mw_se*1.65)
trans_mw_ci_high_pooled_90 <- coef_from_baseline_trans_mw + (trans_mw_se*1.65)

# Plotting
setEPS(width=4,height=4)
postscript("figures/figure2b-final.eps")
par(mar=c(4.5,4.5,1,1), mgp=c(2,.25,0))
plot(coef_from_baseline_gay_mw,c(1.2,0.2),ylim=c(0,2),xlim=c(-.3,.1), tck=0, cex.main=1.4, main="", cex.lab=0.7, cex.axis=0.75, 
     pch=c(16,17), cex=2, xlab="Change in Support from Baseline (Cisgender or Straight)", xaxt="n", yaxt="n", ylab="", lwd=1.5, col=c("grey30","grey60"))
abline(v=0, lty=2)
axis(1,at=seq(-.3,.1,by=.1), cex.axis=0.75, 
     labels=c("-.3","-.2","-.1","0",".1"), tck=0.01, las=1)
axis(2, at=c(0.45,0.15,1.45,1.15),cex.axis=0.75,labels=c("Transgender 
Woman","Lesbian","Transgender 
                                                Man","Gay"), padj=0, las=1, tick=F)
legend("topleft", c("All Respondents, Men Nominees","All Respondents, Women Nominees"), col=c("grey30","grey60"), pch=c(16,17), cex=.7)
segments(gay_mw_ci_low_pooled, c(1.2,0.2), gay_mw_ci_high_pooled, c(1.2,0.2), lwd=1, col=c("grey30","grey60"))
segments(gay_mw_ci_low_pooled_90, c(1.2,0.2), gay_mw_ci_high_pooled_90, c(1.2,0.2), lwd=3, col=c("grey30","grey60"))

points(x=coef_from_baseline_trans_mw, y=c(1.5,0.5), pch=c(16,17), col=c("grey30","grey60"), cex=2)
segments(trans_mw_ci_low_pooled, c(1.5,0.5), trans_mw_ci_high_pooled, c(1.5,0.5), lwd=1, col=c("grey30","grey60"))
segments(trans_mw_ci_low_pooled_90, c(1.5,0.5), trans_mw_ci_high_pooled_90, c(1.5,0.5), lwd=3, col=c("grey30","grey60"))
dev.off()

################################################################
# In-text discussion of coefficients associated with Figure 2b #
################################################################
### In Figure 2b, we illustrate that the magnitude of the negative effect for both transgender
### and gay/lesbian status is larger for women (16.8 and 12.0 percentage points, respectively)
### than men (10.4 and 3.3, respectively)

# Gay penalty for women nominees, coefficient and p-value
round(summary(results.support.just.women)$amce[15,3],3) # coefficient
round(summary(results.support.just.women)$amce[15,6],3) # p-value
# Gay penalty for men nominees, coefficient and p-value
round(summary(results.support.just.men)$amce[15,3],3) # coefficient
round(summary(results.support.just.men)$amce[15,6],3) # p-value

# Transgender penalty for women nominees, coefficient and p-value
round(summary(results.support.just.women)$amce[3,3],3) # coefficient
round(summary(results.support.just.women)$amce[3,6],3) # p-value
# Transgender penalty for men nominees, coefficient and p-value
round(summary(results.support.just.men)$amce[3,3],3) # coefficient
round(summary(results.support.just.men)$amce[3,6],3) # p-value

### this difference is statistically significant for gay/lesbian (p=0.027) but not transgender (p=0.108). 
# Following Clogg, Petkova, and Haritou 1995 (pg.1276)

# Negative effect of gay/lesbian (values come from above models)
b_lesbian <- -0.1196301
b_gay <- -0.0334528
se_lesbian <- 0.026813
se_gay <- 0.028190

z_gay_lesbian <- (b_lesbian - b_gay)/(sqrt((se_lesbian^2)+(se_gay^2)))
round(2*pnorm(z_gay_lesbian*-1, lower.tail = F),3)

# Negative effect of transgender (values come from above models)
b_trans_woman <- -0.1681376
b_trans_man <- -0.1039549
se_trans_woman <- 0.027809
se_trans_man <- 0.028724

z_trans <- (b_trans_woman - b_trans_man)/(sqrt((se_trans_woman^2)+(se_trans_man^2)))
round(2*pnorm(z_trans*-1, lower.tail = F),3)
